clear all; clc;%close all;
%Updated: 24/05/2021

a   =load('agrid2.txt');
b   =load('bgrid2.txt');

def0  =load('gdef1.txt');
gn0   =load('ggn.txt');
gnaut0=load('ggnaut.txt');
dist0 = load('mcdistA2.txt');
index_bor = find(dist0(:,1)==1);
bort0 = dist0(:,1);
at0   = dist0(:,2);
bt0   = dist0(:,3);

at0orig = at0;
bt0orig = bt0;
bort0orig = bort0;
at02 = at0;


na=size(a,1);
nb=size(b,1);

lambda = 0.184;
r      = 0.02;

atrans = ones(nb,1)*a';

bgrid=b;

for i=1:na;
	gnaut(i)=gnaut0(i);   
    for j=1:nb;
        def(i,j)=def0((i-1)*nb + j);   
        gn(i,j)=gn0((i-1)*nb + j);   
        gnstar(i,j)=(1-def(i,j))*gn(i,j)+(def(i,j))*gnaut(i);
    end;
end;

for i=1:na-1;
    for j=1:nb;
        gndif(i,j)=gn(i+1,j)-gn(i,j);  
        gndif(i,j)=100*gndif(i,j);
        if (def(i,j)==0);   %default region
            gndif2(i,j)=2;
        else;
            if gndif(i,j)>0.001;
                gndif2(i,j)=1;
            elseif gndif(i,j)<-0.001;
                gndif2(i,j)=0;
            else;
                gndif2(i,j)=-1;
            end;
        end;
        
    end;
end;

[ndebt,i2]=find(b==0);
gndif3=zeros(na-1,ndebt);
for i=1:na-1;
    for j=1:ndebt;
        gndif(i,j)=gn(i+1,j)-gn(i,j);  
        gndif(i,j)=100*gndif(i,j);
        if (def(i,j)==0);   %default region
            gndif3(i,j)=NaN;
        else;
            gndif3(i,j)=gndif(i,j);
        end;
        
    end;
end;

index_a=find(at02<na);
aat0=at02(index_a);


bt0 = bt0(index_a);
index_b = find(bt0<=ndebt);
bt0 = bt0(index_b);

bbt0=-bgrid(bt0)/(lambda+r);
aat0=aat0(index_b);

gdpss=4.822;
aat02=exp(a(aat0));

bbt02=100*bbt0/gdpss;

i_at0 = aat0;
i_bt0 = bt0;
i_bort0=bort0(index_a);

nsim=rows(i_at0);

i0=0;
i1=0;
i2=0;
im1=0;

for i=1:nsim;
    i_a=i_at0(i);
    i_b=i_bt0(i);
    if (gndif2(i_a,i_b)==2);
        i2=i2+1;
    elseif (gndif2(i_a,i_b)==1);
        i1=i1+1;
    elseif (gndif2(i_a,i_b)==0);
        i0 = i0+1;
    elseif (gndif2(i_a,i_b)==-1);
        im1=im1+1;
    end;
end;

a=exp(a);

b =-bgrid/(lambda+r);  %change of sign
b=100*b/gdpss;
bbt0=100*bbt0/gdpss;
aat0=a(aat0);

bx = b(1:ndebt);
ax = a(1:end-1);

format_chart='epsc';
plot_path=pwd;

nameFig = 'heatmap';figure('name',nameFig);
h=heatmap(bx,ax,gndif3,'GridVisible','off','MissingDataLabel','default region','MissingDataColor',[0.988235294818878 0.882352948188782 0.894117653369904])
h.XDisplayData = flipud(h.XDisplayData);  
h.YDisplayData = flipud(h.YDisplayData);  
 h.FontSize = 14;
bx=flipud(bx);
axx=flipud(ax);
CustomXLabels=string(bx);
CustomXLabels(1:end)= ' ';
h.XDisplayLabels = CustomXLabels;
h.XLabel = 'debt (% GDP^{ss})';
CustomYLabels=string(ax);
CustomYLabels(1:end)= ' ';
h.YDisplayLabels = CustomYLabels;
s = struct(h);
s.XAxis.TickLabelRotation = 0;  
h.XDisplayLabels = CustomXLabels;
ax = axes;
hh=scatter(bbt02,aat02,'filled','k');box off;
axis([bx(1) bx(end)-1 a(1)-0.0075 a(end-1)+0.0075]);
ax.Color = 'none';
ax.XTick = [0, 10, 20 , 30, 40]; 
set(ax, 'FontSize', 14)
ylabel(' y^{T}', 'FontSize', 14);
saveas(gcf, fullfile(plot_path, nameFig),format_chart);